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Cite as : C. Vega, M.Martin-Conde and A.Patrykiejew, Mol.Phys., 104, 3583, (2006). 
Molecular dynamic simulations were performed for ice Ih with a free surface. The simulations were 

'^ carried out at several temperatures and each run lasted more than 7ns. At high temperatures the 

^ ice melts. It is demonstrated that the melting process starts at the surface and propagates to the 

f^ bulk of the ice block. Already at the temperatures below the melting point, we observe a thin 

o 

O liquid layer at the ice surface, but the block of ice remains stable along the run. As soon as the 

,__( temperature reaches the melting point the entire ice block melts. Our results demonstrate that, 

> 

QQ unlike in the case of conventional simulations in the NpT ensemble, overheating of the ice Ih with 

(^ a free surface does not occur. That allows to estimate the melting point of ice at zero pressure. 

in 

^' We apphed the method to the following models of water: SPC/E, TIP4P, TIP4P/Ew, TIP4P/Ice 

^r^ and TIP4P/2005, and found good agreement between the melting temperatures obtained by this 

;v procedure and the values obtained either from free energy calculations or from direct simulations 

KA of the ice/water interface. 



I. INTRODUCTION 

For obvious reasons, liquid water has been the focus of thousands of simulation studies 
since the pioneering work of Barker and Watts^ and Rahman and Stillinger-. However, 
simulation studies of the solid phases of water (ices) have been much more scarce^^. There 
is a number of reasons to study solid and amorphous phases of water by computer simulation 
methods. First of all, the description of the phase diagram represents a major challenge 
for any potential model of water. Secondly, many experimental facts concerning ices and 
amorphous water are not completely understood and computer simulations could help in 
obtaining a molecular view of the process. Just to show a few examples, let us mention the 
nucleation of ice^'^, the dynamics of solid-solid transitions, the possibility of liquid-liquid 
equilibriaP"^, the speed and mechanism of crystal growthP^*^ and the properties of ice at 
the free surfacd^^^^. 

Any researcher performing simulations of water must choose a potential model among the 
many now available^^"^. The most popular were developed in the eighties and are known 
as: TIP3FI231 , TIP4Pl21l25]^ TIPSplSSl, SPCPand SPC/E^s. They were fitted to reproduce the 
properties of liquid water at room temperature and pressure. Before choosing a potential 
model of water it seems quite reasonable to ask about its ability to describe the phase 
diagram of water. Over the last decade, the vapour - liquid coexistence has been determined 
for several water models by using the Gibbs ensemble technique^^nSSl. The critical properties 
of many water models are now known. However, the melting point of these models has been 
studied less often. Haymet et al. provided an estimate of the melting point of SPG, TIP4P 
and SPG/E water models^^'^, by performing simulations of the liquid-solid interface. The 
same approach has been followed recently by Wang et al.'^ and by Nada and FurukawaP'. By 
performing free energy calculations the melting points of the models SPG/E, TIP4P ,TIP5P 
and NvdEp2 have been reportecP'^^"^. It should be stated that estimates of the melting 
point of ice obtained by different authors and different techniques (or even by the same 
authors using different techniques) do not always agree, so it is still of interest to determine 
the melting point by as many different routes and methodologies as possible. Over the 
last few years, we have undertaken a systematic study of the solid phases of water^™^. 
Firstly the phase diagrams for SPG/E and TIP4P models were obtained from free energy 
calculations^SEZI \yy using the methodology proposed by Frenkel and LadcP^, and extended 



to water models by Vega and Monsoii^. Secondly, by using the Hamiltonian Gibbs-Duhem 
integratioEp"^ the melting point of ice Ih was obtained for several models of water "^^. In 
these two last cases we used Monte Carlo methods and "home-made" codes. Then, we 
applied a completely independent methodology, namely molecular dynamics method^, and 
used the GROMACS package^ to perform direct simulations of fluid-solid coexistenc d^^ * ^° l 
It should be emphasized that we obtained excellent agreement between the normal melting 
temperatures T^ ( i.e the melting temperature at the normal pressure oip = Ibar ) estimated 
by different routes. That was the case for the traditional SPC/E, TIP4P and TIP5P, and also 
for the new generation of models proposed just in the last two years, namely TIP4P/Ew^, 
TIP4P/Ice^ and TIP4P/2009^, designed to improve the description of ices and water. Of 
course, each model predicts a different normal melting temperature, and the TIP4P/Ice^ 
model is known to perform best leading to Tm = 270 ± 2K. 

One may wonder why it is so cumbersome to determine the melting point of a water model. 
After all, melting points are determined easily in the lab conditions for many substances 
by simply heating the solid at constant pressure until it melts. Although liquids can be 
supercooled, solid can not be superheated as first stated by BridgmanrpS " It is impossible 
to superheat a crystalline phase with respect to the liquid. ". For this reason there is in 
principle no risk in determining melting points by simply heating the solid until it melts. 
However, in computer simulations the situation is somewhat different. If conventional NpT 
simulations are performed for ice Ih then it does not melt at Tm but at a somewhat higher 
temperatur#^™^. It is worth mentioning that superheating has been also recently found 
experimentally^SEQl fQj^ j^g j^ i^^^ only on a short time scale (about 200ps). The difference 
between the results of NpT simulations and those found in experiments (summarized in 
the Bridgman's statement) is striking. However, there is a fundamental difference between 
experiments and NpT simulations of bulk solids. Whereas in experiments the ice (or a solid 
in general) must have a surface (or rather an interface at which it is in a contact with another 
phase; vapour or liquid), this interface is missing in conventional NpT simulations of bulk 
solids. When no interface is present ice must melt via bulk melting^. 

It is now commonly accepted that melting starts at the surface and already at the temper- 
atures lower than the bulk melting point, solids exhibit a liquid-like layer at the surfac#^^™. 
Different melting scenarios are possible depending on the behaviour of that liquid-like layer 
when the bulk melting temperature is approached from below. When the thickness of the 



liquid-like layer increases and finely diverges at T = Tm, i.e., when the solid is wetted by its 
own melt, one meets the so-called surface induced melting or just surface meltingZ3I21Il3{Z9l_ 
On the other hand when the liquid-like layer retains finite thickness up to T = T^, i.e., when 
the melt does not wet completely the solid surface (partial wetting), one meets the case of 
surface nonmelting^^'^. Both situations have been observed experimentally and found in 
computer simulations (for a review see that of Tartaglino et al.!^. No matter which of the 
two above mentioned mechanisms occurs, the existence of an interface prevents the solid 
from overheating. One may therefore expect that the presence of a free surface at the ice 
sample, would turn the results of computer simulations back to normal, i.e. with melting 
occurring right at the bulk melting temperature. 

The purpose of this paper is two-fold. On one hand we would like to check whether a 
piece of ice with a free surface can or not be superheated in computer simulation. As it 
will be shown shortly, the answer to this question is that when a free surface is present, 
superheating is suppressed in computer simulation, just the same as in experiment. Things 
go back to normal, and Bridgman's statement holds again. Secondly the suppression of 
superheating means that there is another relatively straightforward way to determine the 
normal melting point. As it will be shown , the melting point obtained from the simulation 
of the free surface agrees with the estimates obtained by other routes. The problem of 
estimating melting points of water models seems solved, remaining uncertainties being of 
the order of about 3-4K. 

II. METHODOLOGY 

Let us first describe the procedure used to get the initial configuration. Although the ice 
Ih is hexagonal, it is possible to use an orthorrombic unit cell^^. It was with this orthorrombic 
unit cell that we generated the initial slab of ice. In ice Ih, protons are disordered whereas 
still fulfilling the Bernal-Fowler rules^^'^. We used the algorithm of Buch et alP^ to obtain 
an initial configuration with proton disorder and almost zero dipole moment (less than 0.1 
Debye). This initial configuration contained 1024 molecules of water and fulfilled the Bernal 
Fowler rules. The dimensions of the ice sample were 36 A x 31 A x 29 A . In oreder 
to equilibrate the solid, NpT simulations of bulk ice Ih were performed at zero pressure at 
each temperature of interest. We used the molecular dynamics package Gromacs (version 



3.3p^. The time step was Ifs and the geometry of the water molecules was enforced using 
constrainta^^'^. The Lennard- Jones of the potential (LJ) was truncated at 9.0 A . Ewald 
sums were used to deal with electrostatics. The real part of the coulombic potential was 
truncated at 9.0 A . The Fourier part of the Ewald sums was evaluated by using the Particle 
Mesh Ewald (PME) method of Essmann et al^. The width of the mesh was 1 A and we 
used a fourth other polynomial. The temperature was kept by using a Nose-Hoover^SES 
thermostat with a relaxation time of 2 ps. To keep the pressure constant, a Parrinello- 
Rahman barostatF^'^ was used. The relaxation time of the barostat was of 2 ps. The 
pressure of the barostat was set to zero. We used a Parrinello-Rahman barostat where 
all three sides of the simulation box were allowed to fluctuate independently. The angles 
were kept orthogonal during this NpT run, so that they were not modified with respect to 
the initial configuration. The use of a barostast allowing independent fluctuations of the 
lengths of the simulation box is important. In this way, the solid can attain the equilibrium 
unit cell for the considered model and thermodynamic conditions. It is not a good idea 
to impose the geometry of the unit cell. The system should rather determine it from NpT 
runs. Once the ice is equilibrated at zero pressure we proceed to generate the ice-vacuum 
interface. By convention, we shall assume in this paper that the x axis is perpendicular to 
the ice-vacuum interface. The ice-vacuum configuration was prepared by simply changing 
the box dimension in the x axis from its value ( around 36 A ) to 100 A . As a result, 
the dimensions of the simulation box are now 100 Ax31 Ax 29 A (the y and z are 
approximate since the actual values were obtained from the NpT runs for each model and at 
any particular value of the temperature). The approximate area of the ice- vacuum interface 
was 31 A X 29 A (approximately). This is about 10 molecular diameters in each direction 
parallel to the interface. In this work we have chosen the x axis to be of the 1210 direction. 
In other words, the ice is exposing its prismatic secondary plane to the vacuum (i.e the 1210 
plane). In fig.l we show the initial configuration of a perfect ice as seen from the basal plane. 
More details regarding ther relation of the main planes of ice (basal, primary prismatic and 
secondary prismatic) to the hexagonal unit cell can be found in in figure 1 of the paper by 
Nada and Furukawa^^, in the paper by Carignano et al.'^ and also in the water web site of 
ChaplinP'. Obviously, the properties of the ice- vacuum interface depend on the crystal plane 
facing the solid-vacuum interface. In particular, the melting mechanism (surface melting or 
surface nonmelting) may depend on the structure of the surface plane. On the other hand. 



the melting temperature should not depend on the choice of the plane selected to form the 
ice- vacuum interface. Therefore, the choice of the prismatic secondary plane is a good as any 
other. The main advantage arises from the fact that for the ice-water interfaces it has been 
shown that the 1210 plane exhibits the fastest dynamics^^'^. It is also a convenient face 
to make a visualization of the melting process. Once the ice-vacuum system was prepared, 
we performed relatively long NVT runs (the lengths was between 4ns and 14ns depending 
on the model and thermodynamic conditions). Since we have been using NVT molecular 
dynamics, the dimensions of the simulation box have been fixed, of course, unlike in the 
preceding NpT run. 

III. RESULTS 

Let us first present results for the TIP4P model. In fig.2 the results of the total energy 
along the run are presented. As it is seen the total energy presents two different behaviours 
depending on the temperature. At high temperatures, the total energy of the system in- 
creases continuously and then reaches a plateau. The visualization of particle trajectories 
indicates that the block of ice melts and finally one obtains an slab of liquid in equilibrium 
with its vapour. The melting speed depends on the temperature but as it can be seen, 
close to the model melting point it is possible to melt the ice completely in about 6ns. The 
difference between the initial and the final energy is about 3.5 kJ/mol, indicating that the 
melting enthalpy is roughly of this order of magnitude. This is in reasonable agreement with 
the melting enthalpy of the TIP4P model, estimated in our previous work"^^ for the melting 
of bulk ice to bulk water, which was of about 4kJ/mol (the somewhat smaller value found 
here is due to the presence of the free surface ). Since the block of ice is of a thickness of 
about 36 A , one can approximately state that an increase of energy by about 1 kJ/mol cor- 
responds to a decrease of the thickness of the ice block by about 10 A . Taking into account 
that we have two solid/vacuum interfaces, one at the right and the other at the left side of 
the ice block, it means roughly speaking an increase of the total energy by IkJ/mol, which 
amounts to the formation of a liquid layer of about 5 A (the argument is not elaborated but 
at least provides some orders of magnitude). The behaviour at low temperatures is different. 
At the beginning (first l-2ns), there is an increase of the energy but after that the energy 
remains approximately constant, apart from the thermal fluctuations. The analysis of the 



configurations of tlie TIP4P at T = 228K, sliows tliat tlie increase of energy during the first 
Ins is due to tfie formation of a thin liquid layer at the surface of ice, which may indicate 
the onset of surface melting, mentioned already in the Introduction, and first proposed by 
Farada}^. 

We just recall that the term "surface melting" usually represents the formation of a quasi- 
liquid layer (qll) on the surface of a solid at temperatures still below the melting point. The 
properties of this quasi-liquid layer are similar to, but not identical with, those of a bulk 
fluid under the same conditions. The thickness of the layer depends on the particular 
plane forming the solid-fluid interface and on the proximity to the melting point. It may 
diverge to infinity or stay finite as the temperature approaches T^ from below. The surface 
melting of ice, has been found both in experiment (see^^'^ and references therein ) and 
in computer computer simulation for several potential models of water^^'^. Some theories 
are also available to explain its origiiF^Ml One should note that the surface melting is an 
equilibrium phenomenon, leading to the lowering of the system free energy by replacing 
the ice-vacuum interface by the ice-liquid* and liquid*-vapour interfaces. By liquid* we 
mean a quasi-liquid layer of microscopic thickness. The existence of a quasi-liquid layer 
at temperatures below the melting point is not exclusive of ice, but it does also appear in 
simpler liquids such as metals^^'^ or Lennard- Jones models^^'^. The study of the quasi- 
liquid layer thickness by computer simulation is of interest by its own. This is so because 
different experimental techniques provide quite different values of its thickness^^'^ and it 
is of interest to understand why this is the case. We expect to study that in future work. 
For the purpose of this work it is just enough to state, that at low temperatures a qll 
appears at the ice surface, provoking an increase of energy in the first l-2ns of the run, 
but then the energy remains constant (with the normal thermal fluctuations) after this 
initial period. In that respect no signal of melting of the block of ice has been observed 
at the temperatures below Tm- By repeating the simulation at several temperatures it is 
possible to determine the lowest temperature at which the block of ice melts T+, and the 
highest temperature at which it does not melt T_. By taking the average of these two 
temperatures we obtain what we call T^ = 1/2(T_|_ + r_). For the length of the runs used 
(about 4-lOns) the typical difference between T+ and T_. is about 3-4K. It would be possible 
to reduce this temperature window by performing longer runs (of the order of hundred ns 
rather than of 10ns). However, taking into account that the current uncertainty in the 



estimations of the water melting point for different models obtained from the free energy 
calculations is just about 4K, the accuracy obtained by the runs of the length presented 
here seems to be quite sufficient. The crucial question now is, what T^ actually is? How 
does Ts compare with the melting point T^,? For the case of TIP4P, the melting point 
has been determined by several groups. From the free energy calculations we obtained^, 
Tm = 232K. Starting from the melting point of the SPC/E model obtained by free energy 
calculations ( T^ = 215K ) and using Hamiltonian Gibbs-Duhem integration we obtained 
again T^ = 232K for the TIP4P modeW Using free energy calculations, Koyama et alP^ 
have obtained T^ = 229K. From direct coexistence between the fluid and the solid phase we 
have obtained Tm = 228i^d^. From direct coexistence between the fluid and the solid Wang 
et al."^ have obtained Tm = 229K. From the results presented in four different papers, each 
using different techniques and methods, it seems clear that the melting point of the TIP4P 
model is about T = 230 (3) -ft'. This is also consistent with a corresponding states rule found 
by us recently. We have founcpSl that for TIP4P model the ratio of the melting point to the 
critical point temperatures is about 0.394. This estimate of the melting point temperature 
is consistent with the corresponding states rule and with the well known value of the critical 
temperature of the TIP4P modeP^'^^'^. The fact that there are some small discrepancies 
between different estimates is not surprising. They are due to the fact that long range 
forces are treated in a different way in different papers (truncation of the potential versus 
Ewald sums), different ice configurations (recall the issue of the proton disorder), different 
cutoff of the potential, as well as possible finite-size effects. Therefore, it can be stated with 
confidence that the melting point of the TIP4P is about T = 230 (3) fC. The value of Tg 
found in this work is Tg = 229(2)^. The error in Tg is mainly determined by the grid of 
temperature used, being the difference between T+ and T_ of about 4K, so that the error 
in Tg is of about 2K. The obvious conclusion is that for the TIP4P model, and within the 
uncertainty of different estimates, Tg is identical to Tm- In other words, overheating has 
been completely suppressed by the presence of a free surface. 

In order to test the hypothesis, stating that T, and Tm are identical, we have per- 
formed runs for other different models of water, namely the recently proposed TIP4P/Ice^ 
, TIP4P/200dS31 , TIP4P/EW5I] and the classical SPC/E model. Resuhs of the runs for 
these models are presented in Figs. 3-6 respectively. In Table I, the values of Tg for different 
models are presented and compared with the values of Tm obtained from a direct study of 
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a fluid-solid coexistence, as well as from the Hamiltonian Gibbs-Duhem integration. As it 
can be seen, T, is identical to T^, within the error bar. This result seems to be impor- 
tant, since it shows that the similarity between Tg and T^ is not a particular feature of the 
TIP4P model, but applies to any other model of water. The conclusion is that when a free 
surface is available overheating of the solid is suppressed, and the ice melts at the melting 
temperature, provided the length of the runs are sufficient (of the order of 10ns or larger), 
and the size of the ice block is sufficiently large, not smaller than about 3-4nm. To check 
whether this conclusion depends on the system size and/or of the face exposed to vacuum 
we have performed a run for the TIP4P/Ice model at a temperature of T=274K (higher 
than the melting point). In this particular case we used 1536 molecules and the plane at 
the ice-vacuum interface was the primary prismatic plane. The area of the interface was of 
about 27 A X 30 A . The width of the ice was of 63 A approximately. Again, complete 
melting of the sample was observed, indicating that the physics was not changed neither by 
a larger system nor by a different plane exposed. However, in this case it took about 25ns 
to melt the ice completely due to a larger size of the system and to slower dynamics. 

A minor comment is in order here. The melting point obtained in our previous work T^ 
corresponds to the melting point at the normal pressure p = Ibar. However, Tg was obtained 
in this work from simulations at zero pressure. Obviously, the melting temperature at zero 
pressure and at the pressure of one bar differ somewhat. In fact, for real water the difference 
is about O.OIK (i.e 273. 16K for essentially zero pressure at the triple point versus 273. 15K for 
p = Ibar). For water models, the difference between the melting pressure at zero pressure 
and at p = Ibar is expected to be of the order of l/{dp/dT). By using the values of dp/dT 
reported in our previous workF^ it can be seen that also for water models the difference 
between the melting temperatures at normal pressure and zero pressure should be of the 
order of about O.OIK (i.e quite small). In view of that, it it reasonable to identify Tm with 
T,. Still another comment is required with respect to our simulations of the ice- vacuum 
interface. In principle, when ice is exposed to vacuum at temperatures below the melting 
point there should be a solid-vapour equilibria. The vapour pressure of water at the triple 
point is quite low, namely 6 x 10~^ bar. The vapor pressure of TIP4P models at the triple 
point is fifty times smaller^, about 1.3 x 10~^ bar. Although in principle one could determine 
the vapour pressure of ice by performing NVT runs of ice in contact with vacuum, this is 
not feasible in practice. In fact, for the lengths of the runs used in this paper we have never 



observed the sublimation of a single ice molecule into the vacuum (this of course should 
occur in very long runs, but since it is a rare event it was not observed in the runs of this 
work which lasted about 10ns). For this reason we prefer to say that T, was obtained at zero 
pressure, but it could be more appropriate to say that it was obtained at the sublimation 
pressure of the ice (which must be quite low anyway) . 

It is interesting to see the mechanism of the melting process. We shall present results for 
the TIP4P/Ice at T = 2'IQK (a temperature higher than the melting point of the model). 
In fig. 7, a snapshot of the system after Ins is shown. In fig.8, the final snapshot of the 
system after 5ns of run is shown. As shown in fig.8, the ice has melted completely after the 
run of 5ns. From the results given in Fig. 7 it is also clear that the melting process starts 
at the surface and then propagates to the interior or the block of ice. We have analyzed 
the movies of all the runs and always found this to be the case. We have never observed 
the formation of a droplet of water in the interior part of the ice. We have always observed 
that the melting started at the surface and propagated to the interior part of the block 
of ice. One of the first methodologies proposed to estimate the melting point in computer 
simulations was the direct coexistence metho(J ^^ * ^° * ^^ l In this method, the solid is put into 
contact with the liquid at a certain temperature and pressure. If the temperature is higher 
than the melting temperature, the solid part of the sample will melt. If the temperature 
is lower than the melting temperature then the liquid will freeze. Regarding the ice- water 
interface, it has been investigated during the last fifteen years and estimates of the melting 
points for several water models have been given.^^^^ ^^^*^^*^^ i Recently, we have found that 
melting temperatures estimated from the direct coexistence method are in agreement with 
those obtained from free energy calculations. The presence of a "nucleus" of ice in the 
sample, and of a "nucleus" of liquid in the sample prevents the possibility of met ast ability. 
In other words, the presence of the two phases in the simulation sample from the very 
beginning avoids the possibility of metastability. The metastability is due to the fact that 
the formation of a nucleus of a new phase in the interior of another phase is an activated 
process. Although this is pretty clear, it was not so obvious what happened when the free 
surface of the solid was heated while exposed to vacuum. This work shows that also in 
this case superheating is suppressed. To prove further that a liquid layer is already present 
in the sample at temperatures well below the melting point, we present in fig. 9 the final 
configuration (after a 8ns run) obtained for the TIP4P/Ice at a temperature well below 
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the melting point of the model (T = 2QAK). As it can be seen, a quasi liquid layer is 
already present in the system. In fig. 10 the final configuration (after a 9ns run) obtained for 
the TIP4P/Ice at a temperature just below the melting temperature T = 270-/^ is shown. 
It can be observed that the thickness of the quasi-liquid layer is larger now that at the 
lower temperature. The question whether the thickness of the liquid layer diverges when 
the melting point is approached (total wetting) or remains finite (partial wetting) should 
be studied in more detail. In any case figures 9 and 10 illustrate clearly the point that a 
quasi-liquid layer is already present in ice at temperatures below the melting point so that 
superheating is suppressed. 

IV. CONCLUSIONS 

In this work, we have performed NVT Molecular Dynamics simulations of ice Ih with a 
free surface. The ice was first equilibrated by performing NpT simulations at zero pressure. 
Simulations were then performed at several temperatures and runs lasted between 5 and 
10ns. At temperatures below the melting point, the ice surface develops a thin liquid like 
layer (quasi liquid layer) of microscopic thickness. However, at temperatures above the 
melting point the ice melt, and does not exhibit any trace of overheating. In all cases, we 
have observed that the melting process starts at the surface and propagates to the interior 
of the material. We have estimated the melting point of ice Ih for several water models, 
namely , TIP4P, TIP4P/Ice, TIP4P/2005, TIP4P/Ew, and SPC/E. Since the mefting points 
of these water models had been evaluated previously (for ice Ih) by free energy calculations 
and by fiuid-solid interface simulations, it is of interest to compare to the values obtained in 
this work. Excellent agreement has been found (see Table I) between previous estimates of 
Tm and those obtained in this work from simulations of the free surface. Therefore, although 
it is possible to overheat ice in the usual BULK NpT simulatiorP^'^, the overheating does 
not seem to be possible as soon as the ice has a free surface. Although this phenomenon has 
been studied before for spherical fiuidsSDlQQl^ to the best of our knowledge this is the first 
time that the issue is addressed for such a complex fiuid as water. As a bonus, one has a 
remarkably simple method to estimate the melting temperature of a water model: equilibrate 
first the material in an NpT run, and perform afterwards a long NVT run (of about 10ns 
or more). This is a long run but still within the capacity of current computers. Taking 
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into account that estimating melting temperatures by free energy calculations is somewhat 
involved (although feasible) the existence of an alternative and simpler method is welcome. 
One should note that the method is probably not limited to water, but can also be used for 
other complex molecules. The method can be applied only when the following two conditions 
are met. Firstly, the melting of the system at temperatures above the melting point should 
be fast enough to be studied within a reasonable time by molecular simulations. Secondly, 
the liquid layer thickness at the temperatures below, but not too far from the melting point, 
should not be too large. Otherwise, the finite size of the sample may lead to the observation 
of complete melting already below the melting point. It should be pointed out that the 
procedure allows to estimate the zero pressure melting temperature, which should be very 
close to the normal melting point, but does not allow to estimate the fluid-solid coexistence 
at higher pressures. 

Now, back to the classroom you may explain that the fluid-solid transition (in the ab- 
sence of impurities leading to heterogeneous nucleation) occurs by the following mechanism: 
nucleation (of an embryo of the solid phase) and growth. Liquid water can be supercooled 
because the formation of an embryo is an activated process and that requires a certain 
amount of time to occur. You may now explain that melting occurs by the same mech- 
anism: nucleation (of am embryo of the liquid phase) and growth. The point is that for 
ice with a free surface the embryo of the liquid phase is already there, at the surface of 
the solid. For this reason, the activation energy of forming the liquid embryo is zero and 
melting is just the growth of a new phase (i.e the propagation of the liquid layer from the 
surface to the interior of the solid). As stated by Bridgmann in his classical 1912 paper^ " 
It is impossible to superheat a crystalline phase with respect to the liquid. No good reason 
for this has ever been given, but no exception has ever been found, and it is coming to 
be regarded as a law of nature " . We agree except for the absence of explanation for this 
behaviour. As stated by FrenkeP^'^^ and suggested by Tammann|i^ " It is well known that 
under ordinary conditions an overheating of a crystal, similar to the overheating of a liquid, 
is impossible. This peculiarity is connected with the fact that the melting of a crystal , 
which is kept at a homogeneous temperature, always begins on its free surface. The role of 
the latter must, accordingly, consist in lowering the activation energy, which is necessary for 
the formation of a flat embryo of the liquid phase, i.e. of a thin liquid film down to zero." 
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TABLE I. Normal melting points Tm of different water models as obtained from different method- 
ologies. The column labeled as free energy correspond to the melting point either as obtained from 
free energy calculationP^EII (SPC/E and TIP4P) or from Hamiltonian Gibbs Duhem integratiorP^. 
(TIP4P/Ice, TIP4P/2005, TIP4P/Ew). The results of the column labeled as Sohd-Liquid interface 
have been obtained from direct coexistence of the solid with the liquid^. The last column is the 
value of Tg (which is just the average of the lowest temperature at which ice melts and the highest 
at which it does not melts) obtained from simulations of ice Ih with a free interface. The similarity 
of Ts and Tm shows the absence of superheating for ice with a free surface. 

Model Tm/K (Free energy) Tm/K (Solid-Liquid interface) Tg/K (this work) 

229(2) 230(2) 

268(2) 271(1) 

249(2) 249(3) 

242(2) 243(2) 

213(2) 217(2) 



TIP4P 


232(4) 


TIP4P/Ice 


272(6) 


TIP4P/2005 


252(6) 


TIP4P/EW 


245.5(6) 


SPC/E 


215(4) 
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FIG. 1. Initial configuration of the ice Ih witli 1024 molecules. The x axis goes from left to right 
so that the ice-vacuum interface is located at the right and left sides. The plane that can be seen 
is the basal plane. The plane exposed to the vacuum is the secondary prismatic plane. 
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FIG. 2. Total energy as a function of time obtained at several temperatures by performing MD 
simulations for a block of ice Ih with a free surface. Results presented here correspond to the 
TIP4P model. 
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FIG. 3. Total energy as a function of time obtained at several temperatures by performing MD 
simulations for a block of ice Ih with a free surface. Results presented here correspond to the 
TIP4P/Ice model. 
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FIG. 4. Total energy as a function of time obtained at several temperatures by performing MD 
simulations for a block of ice Ih with a free surface. Results presented here correspond to the 
TIP4P/2005 model. 
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FIG. 5. Total energy as a function of time obtained at several temperatures by performing MD 
simulations for a block of ice J/j with a free surface. Results presented here correspond to the 
TIPAP/Ew model. 
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FIG. 6. Total energy as a function of time obtained at several temperatures by performing MD 
simulations for a block of ice Ih with a free surface. Results presented here correspond to the 
SPC/E model. 
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FIG. 7. Instantaneous configuration of tlie TIP4P/Ice system at T = 276K after Ins run. 
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FIG. 8. Instantaneous configuration of the TIP4P/Ice system at T = 276K at tlie end of the run 
after 5ns. 
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FIG. 9. Instantaneous configuration of the TIP4P/Ice system at T = 264i^ at the end of a 8ns 
run. Although the temperature is well below the melting point of the model, a quasi-liquid layer 
is clearly present in the ice-vacuum interface. 
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FIG. 10. Instantaneous configuration of tlie TIP4P/Ice system at T = 270K at tlie end of a 9ns 
run. The temperature is just below the melting temperature of the model. A quasi-liquid layer is 
clearly present in the ice-vacuum interface. 
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